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13  ABSTRACT 


The  primary  objective  of  the  current  effort  was  the  development  and  solution  of  a  nonlinear  analytical 
longitudinal  Instability  model,  which  would  allow  all  of  the  various  governing  phenomena  to  be  accounted 
for  In  a  coupled  manner.  The  r.\o  primary  elements  of  the  current  instability  analysis  are  a  method  of 
characteristics  solution  of  the  tvo  phase  flow  in  the  combustion  chamber  of  the  motor,  and  a  coupled  calcu¬ 
lation  of  a  transient  burning  rate.  The  transient  burning  rate  analysis  presented,  herein.  Is  a  unique  and 
Interesting  development.  It  Is  based  on  an  extension  of  the  most  popular,  linear,  harmonic  combustion 
response  model.  The  current  method  allows  the  calculation  of  propellant  burning  response  to  a  pressure 
disturbance  of  arbitrary  waveform,  for  all  time,  including  the  period  immediately  following  the  initiation 
of  the  disturbance.  The  analysis  also  includes  a  model  for  velocity  coupled  response.  Therefore,  for 
the  first  time,  the  nonlinear  effects  of  velocity  coupling  on  the  growth  of  pressure  waves  In  a  combustion 
chamber  can  be  computed. 

The  instability  solution,  itself,  begins  with  the  calculation  of  the  steady  state  two-phase  flow  In  the  motor. 
The  flow  In  the  combustion  chamber  is  calculated  by  numerically  integrating  tne  equations  of  motion.  The 
nozzle  flow  is  found  using  the  constant  fractional  lag  approximation.  The  steady  state  conditions  are  then 
perturbed  and  the  subsequent  wave  motion  in  the  motor  Is  calculated  numerically,  using  the  method  of 
characteristics.  The  nature  of  the  engine  response  is  dependent  upon  the  Interaction  the  various  gain  and 
loss  mechanisms  in  the  crglno,  which  are,  in  turn,  a  function  of  the  propellant  burning  response,  the  size 
and  amount  of  particulate  matter  present,  the  magnitude  and  shape  of  the  initial  disturbance  and  the 
geometrical  configuration  of  the  motor. 

The  instability  model  Is  currently  subject  to  the  following  limitations.  Only  motors  with  cylindrlcally 
perforated  grain  w«re  considered.  The  aasdynamic  flew  was  assumed  to  be  one-dimensional  and  the 
particles  in  the  qas  stream  were  taken  to  be  of  uniform  size  and  inert.  The  nozzle  flow  Is  assumed  to  be 
quasi-steady. 

A  scries  of  Instability  solutions  have  been  calculated,  — *  aln  some  of  the  main  parameters  such  as 
particle  size,  burning  rate  constants,  and  Initial  dlst  .  •  .,e  waveform  and  magnitude  have  been  varied, 
in  an  attempt  to  qualitatively  assess  the  behavior  and  Aridity  of  the  present  model.  From  all  appearances, 
the  behavior  of  the  model  13  quite  realistic  and  ltmlted  comparisons  with  data  have  been  quite  encouraging. 
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FOREWORD 


The  present  report  is  part  of  a  two  volume  set  which  describes 
a  nonlinear  solid  rocket  motor  instability  analysis  and  computer  pro¬ 
gram.  Volume  I  contains  the  analytical  basis  for  the  computer  program 
and  a  discussion  of  the  results  obtained  to  date:  Volume  II  of  the  set 
describes  the  computer  program  and  serves  as  a  user's  manual. 

This  investigation  is  entitled  NUMERICAL  ANALYSIS  OF 
NONLINEAR  LONGITUDINAL  COMBUSTION  INSTABILITY  IN  METALLZED 
PROPELLANT  SOLID  ROCKET  MOTORS.  The  two  volumes  are  additionally 
subtitled  as  follows: 

Volume  I  -  Analysis  and  Results 

Volume  II  -  Computer  Program  User's  Manual 

This  investigation  was  sponsored  by  the 

Air  Force  Rocket  Propulsion  Laboratory 

Director  of  Laboratories 

Edwards,  California  93523 

Air  Force  Systems  Command,  United  States  Air  Force 
under  contract  number  F04611-71-C-0060  with  Robert  J.  Schoner  as 
technical  monitor.  Jay  N.  Levine  of  Ultrasystems  (formerly  Dynamic 
Science)  was  program  manager. 

This  technical  report  has  been  reviewed  and  is  approved. 

Paul  J.  Daily,  Lt.  Col.  USAF 
Chief,  Technology  Division 
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NOMENCLATURE 


burning  rate  parameter,  Eq.  (4-14) 

admittance  function 

nozzle  admittance  function 

gas  only,  sound  speed 

sound  speed  based  on  Pp  and  Tp 

burning  rate  parameter,  Eq.  (4-32) 
also,  fractional  lag  parameter,  Eq.  (5-6) 

burning  rate  parameter  for  velocity  coupling 

ratio  of  solid  to  gas  specific  heats,  C/C 

s  P 

constant  in  steady  state  burning  rate,  Eq.  (3-5) 

particle  drag  coefficient 

erosive  burning  constant,  Eq.  (3-5) 

see  el 

specific  heat  of  gas  at  constant  pressure 
specific  heat  of  solid  particles 
specific  heat  of  gas  at  constant  volume 

defined  in  Eq,  (7-34) 
port  diameter 

normalized  surface  activation  energy,  Ew/RnT 
also,  integral  defined  by  Eq.  (8-16)  *  ' 

activation  energy  of  surface  reaction 

internal  energy 

particle-gas  interaction  force  per  unit  volume,  Eq.  (3-8) 
frequency 

also,  as  defined  by  Eq.  (8-14) 
defined  by  Eq.  (7-8) 
defined  by  Eq.  (6-17) 
defined  by  Eq.  (4-20) 
enthalpy 

defined  by  Eq.  (8-12) 
fractional  lag  constant,  Eq.  (5-3) 
chamber  fractional  lag  constant,  Eq.(7~44) 
thermal  conductivity,  also  complex  wave  number 


ill 


k  -  thermal  conductivity  of  the  solid  particles 

s 

L  -  length  of  the  grain,  also  fractional  lag  constant,  Eq.  (5-4) 

L'  -  chamber  fractional  lag  constant,  Eq.  (7-44) 

1  -  perimeter  of  the  grain 

M  -  Mach  number,  also  number  cf  points  on  initial  line 

-  Mach  number  at  burning  surface 

iti  -  particle  mass,  also  surface  mass  flux 

Nu  -  Nusselt  number 

n  -  pressure  exponent  in  steady  state  burning  rate 

r>v  -  constant  :,r.  velocity  coupled  analysis,  Eq.  (4-75) 

n  -  exponent  oressu.ro  dependence  of  surface  reaction  rate 

w 

P  -  pressure 

-  reference  pressure  ir  stv-idy  state  burning  rate 

Pp  -  chamber  pressure 

Pr  -  Prandtl  number 

p  -  pressure,  also  need  for  Laplace  transform  variable 

p^  -  defined  by  Equation  (0-10) 

Q{  -  heat  release  per  unit  mass 

Q  -  particle-gas  heat  transfer  rate  per  unit  volume,  Eq.  (3-14) 

C)  -  heat  of  reaction  for  processes  at  burning  surface 

w 

q  -  see  i 

R  -  gas  constant,  also  normalized  throat  radius  of  curvature 

Rq  -  universal  gas  constant 

R^  -  response  function,  Eq.  (4-35) 

Re  -  Reynolds  number  based  on  particle  diameter  and  particle-gas 

relative  velocity 

RHS  -  right  hand  side  of  a  characteristics  compatibility  relation 
r  -  linear  burning  rate 

-  area  of  burning  surface 

s  -  dimensionless  Laplace  transform  variable,  = 

T  .  -  temperature 

Tp  -  adiabatic  flame  temperature 

t  -  time 

t  -  defined  by  Eq.  (7-33) 

C 

u  -  axial  velocity 

ut  -  threshold  velocity 

W  -  defined  by  Equation  (4-28) 

iv 


:1'  e2 


Ci  ,  c 


1'  2 


reaction  rate  divided  by  gas  density 

axial  distance 

growth  constant 

particle  damping  constant 

defined  in  Eq.  (8-35) 

particle  to  gas  weight  flow  ratio 


ratio  of  particle  to  gas  mass  burning  rates,  w  /w 
ratio  of  gas  specific  heats,  Cp/Cv 
a  small  increment  in  time 


equal  to  t  6 


convergence  criteria  for  characteristics  calculations,  also  used 

in  velocity  coupling  analysis  (Eq.  4-/‘v; 

thermal  diffusivity  of  the  propellant 

defined  by  Eq.  (4-27) 

complex  function  of  frequency,  Eq.  (4-8) 

viscosity 

equals  rxA 

s 

density 

density  based  on  Pp  and  Tf 
density  of  the  metal  oxide  particles 
density  of  the  solid  propellant 
particle  radius 
defined  by  Eq.  (4-62) 


nondimensional  time,  rat/4*s,  also  used  in  Section  2  to 
denote  period  of  oscillation 


characteristic  relaxation  time  for  particle  velocity,  Eq.  (3-1) 
characteristic  relaxation  time  for  particle  temperature,  Eq.  (3-1) 
defined  by  Eq.  (5-10) 
phase  angle 

nondimensional  frequency,  Eq.  (4-9) 


mass  burning  rate,  per  unit  length,  per  unit  cross-sectional 
area,  Eq.  (3-4);  also  occasionally  used  for  angular  frequency 
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e  -  end  of  the  propellant  grain 

f  -  flame 

g  -  gas 

i  -  for  the  £th  mode  of  oscillation 

p  -  particle 

o  -  initial  or  stagnation  value 

t  -  at  the  nozzle  throat 

w  -  at  the  burning  surface  of  the  propellant 


Superscripts 


(  )*  -  in  Sections  3,  5  and  6  only,  denotes  a  dimensional  variable 

(  )'  -  denotes  fluctuation 

(  )  -  in  Section  4  denotes  steady  state  variable,  in  Section  5 

denotes  an  "equivalent"  gas  value,  in  Section  7  denotes 
an  average  quantity 

(  )  -  pertaining  to  a  right  running  characteristic 

(  )“  -  pertaining  to  a  left  running  characteristic 

(  )^  -  real  part  of 
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1.  INTRODUCTION 


This  report  describes  and  serves  as  a  user's  manual  for  a  new 
nonlinear  longitudinal  solid  rocket  motor  instability  program.  This 
program  is  capable  of  treating  metallized  propellants  and  considers 
both  pressure  and  velocity  coupled  transient  burning  response. 

The  analysis,  and  assumptions,  upon  which  this  program  is 
based,  are  discussed  in  Volume  1  of  this  report. 

This  volume  contains  a  description  of  all  of  the  subroutines 
and  program  functions  used  in  the  instability  program.  A  dictionary 
of  all  of  the  variables  contained  in  blank  or  labeled  common  is  also 
included.  Descriptions  of  the  input  required  to  operate  the  program, 
and  the  output  generated,  are  presented;  followed  by  a  sample  case 
which  serves  to  illustrate  the  text,  as  well  as  providing  a  method 
for  checking  out  the  program. 


2.  PROGRAM  DESCRIPTION 


Functionally,  the  Nonlinear  Longitudinal  Solid  Rocket  Motor  Instability 
Program  can  be  considered  to  consist  of  two  distinct  sets  of  subroutines: 
those  required  to  solve  for  the  steady  state  flow  in  the  chamber  and  nozzle; 
and  those  which  solve  the  actual  instability  problem.  All  of  the  subroutines 
and  program  functions  which,  together,  comprise  the  Instability  Program,  are 
briefly  discussed,  in  alphabetical  order,  below: 

SUBROUTINE  AUXSUB 

Part  of  the  nest  of  routines  for  solving  the  steady  state  problem.  This 
routine  evaluates  the  derivatives  in  the  differential  equations 


by  calling  subroutine  GFUNC  to  evaluate  the 


SUBROUTINE  BVP  (NEWTL  LEGS.  AUXSUB.  FDET,  IREAD) 


Subroutine  BVP  and  its  associated  subroutines  AUXSUB,  FCALC,  FDET, 
GFUNC,  LEGS,  NEWH,  and  RKAM,  constitute  a  general  package  for  the  solution 
of  two  point  boundary  value  problems.  This  package  is  part  of  the  Dynamic 
Science  library  of  computer  programs  and  detailed  information  on  any  of  these 
subroutines  is  available  on  request. 

Subroutine  BVP  is  the  overall  control  routine  for  this  package.  In  the  cur¬ 
rent  usage,  nominal  values  for  all  of  the  quantities  usually  input  to  BVP  have  been 
set.  This  allows  the  user  to  employ  the  program  without  becoming  familiar  with 
the  whole  complex  BVP  package.  In  order  to  provide  some  measure  of  flexibility, 
however,  the  option  to  modify  several  of  the  input  parameters  has  been  retained, 
through  the  IREAD  option.  (See  the  discussion  in  the  description  of  input  section). 
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FUNCTION  CNU3ELT  (RE.  PR,  SM) 

This  function  computes  the  Nusselt  number  as  a  function  of  Reynolds  num¬ 
ber  (RE)  and  Prandte  number  (PR).  SM  in  the  calling  sequence  is  not  currently 
used. 


SUBROUTINE  D0PL0T 

This  routine  calls  subroutine  MPL0TS  to  plot  the  pressure,  velocity,  par¬ 
ticle  velocity  and  burning  rate  perturbations  versus  axial  distance. 


FUNCTION  DRAG  (RE) 

This  function  computes  the  particle  drag  coefficient  as  a  function  of 
Reynolds  number  (RE). 


SUBROUTINE  FCALC  (AUXSUB,  FDET) 


Part  of  the  nest  of  subroutines  for  solving  the  steady  state  two  point  boundary 
value  problem.  This  subroutine  controls  the  integration  of  the  differential  equations 
and  the  print  out  of  the  solution.  (Controls  the  printing  of  the  conserved  quantities 
und  their  derivatives.  The  printout  for  the  actual  flow  variables  is  controlled  in 
Subroutine  FLOVARS) .  The  printout  is  controlled  through  the  ITRIG  and  L'.CRxG  par¬ 
ameters,  as  discussed  in  the  description  of  input. 


SUBROUTINE  FDET 


—  "  I 

< 

! 

Part  of  the  nest  of  subroutines  for  solving  the  steady  state  two  point  I 

boundary  value  problem.  This  subroutine  calculates  the  difference  between  I 

: 

the  mach  number  at  the  end  of  the  grain  as  calculated  by  integrating  the  equations  | 

of  motion  and  the  mach  number  calculated  from  the  fractional  lag  nozzle  solution.  j 

and  then  prints  out  the  result. 


SUBROUTINE  FLOVARS 


This  subroutine  solves  for  the  flow  variables  themselves  (U,  T,  P,  U  .  T  . 

*  P#  ***p  • 

Op  from  the  conserved  quantities,  as  required  by  the  steady  state  solution  pack¬ 
age. 


The  flow  variables  are  also  stored  according  to  the  mesh  size  used  In  the 
characteristics  solution.  (The  variables  are  stored  every  NTRIG-th  step).  These 
variables  are  then  printed  out,  and  some  are  later  modified  in  Subroutine  LINEX. 

The  steady  state  solution  is  achieved  by  integrating  the  chamber  pressure  until 
the  mach  number  at  the  end  of  the  grain  matches  that  calculated  in  Subroutine 
STEDYST.  The  flow  variables  at  the  end  of  each  iteration  are  printed  out. 


FUNCTION  FORCE  (R,RP.  U,  UP,  T,  TP,  PF) 

This  function  computes  the  drag  force  between  the  particles  and  gas  as  a 
function  of  local  gas  and  particle  properties. 


•■'•v-yK'-w  ™  V'.T 


SUBROUTINE  FPT  (I,  L  KP) 

This  subroutine  solves  the  field  point  unit  process  which  makes  up  the 
bulk  of  the  method  of  characteristics  solution.  The  arguments  I,  J,  KP  cor¬ 
respond  to  points  1,  2  and  3,  respectively,  in  the  analytical  description  of  the 
field  point  solution. 

The  subroutine  faithfully  follows  the  analytical  description  and  except  for 
the  following  points,  does  not  require  further  discussion. 

The  total  burning  rate  at  point  3  (W3)  is  found  by  extrapolating  from  the  last 
time  a  transient  burning  rate  calculation  was  made  (TIMTB(NK)),  using  the  burn  rate 
derivatives  (DWD'7’).  Linear  interpolation  is  used  to  obtain  the  burning  rate  at 
X  locations  between  those  at  which  the  bum  rate  is  found  directly. 

In  the  solution  for  the  particle  density  at  point  3  (RHOP3)  it  was  convenient 
to  drop  the  area  terms  at  this  time.  In  the  future,  if  variable  area  ports  are 
considered,  this  portion  of  the  program  will  have  to  be  modified. 


SUBROUTINE  FRAKLAG  (E.GBAR) 


This  subroutine  calculates  a  fractional  lag  solution  for  the  nozzle.  In  do¬ 
ing  so  it  establishes  compatible  values  of  K  (the  fractional  lag  parameter  and 
GEAR  (the  isentropic  exponent  of  the  equivalent  ideal  gas).  The  parameter  E, 
which  relates  the  actual  gas  mach  number  to  the  equivalent  gas  mach  number,  is 
also  calculated.  At  the  conclusion  of  the  subroutine  the  results  of  the  fractional 
lag  solution  are  printed  out. 
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SUBROUTINE  GFUNCT  (XXI,  F,  G,  PF) 


Ihis  subroutine  is  part  of  the  set  of  steady  state  routines  and  it  evaluates 
the  derivatives  in  the  differential  equations  (see  subroutine  AUXSUB).  The 
variable  XXI  is  the  local  value  of  X  at  which  the  derivatives  are  to  be  evaluated. 


FUNCTION  HEAT  (U.  UP,  R  RP.  T,  TP,  PF) 

This  program  function  calculates  the  heat  transfer  between  the  particles 
and  gas  as  a  function  of  gas  and  particle  properties! 


SUBROUTINE  INPUT 


Ihis  subroutine  performs  the  following  functions: 

1 .  Zeroes  COMMON/CHAR/ 

2.  Sets  nominal  values  for  some  of  the  $DAIA  input  variables  and  some  other 
constants. 

3.  Reads  Title  Card. 

4 .  Reads  $  DATA  Input 

5.  Calculates  additional  initial  values 

6.  Writes  out  the  $  DATA  input  quantities 

7.  Changes  the  units  of  various  quantities  from  those  input  to  those  used  in¬ 
ternally,  as  follows: 
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A  cycle  of  the  characteristics  solution  is  then  computed  in  two  sweeps 
by  repetitively  calling  subroutine  TPT.  The  characteristics  c.re  tested  for  cros¬ 
sing.  If  crossing  occurs,  the  solution  is  terminated.  The  characteristics  cycle 
is  then  completed  by  calling  the  left  hand  and  right  hand  boundary  point  routines, 
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LHB  and  RHB,  respectively. 

At  the  end  of  each  of  the  above  cycles  the  mesh  point  having  the  smallest 
value  of  time  (not  counting  the  right  hand  boundary  point)  is  identified.  Sub¬ 
routine  INTERP  is  then  called  to  interpolate  the  results  back  into  a  rectilinear 
mesh,  at  the  original  x  locations  and  time  =  t  min. 

Every  NBCALC  characteristics  cycles,  a  new  transient  burning  rate  calcula¬ 
tion  is  obtained.  After  calling  subroutine  TRBRNA  to  update  counters,  subroutine 
TRBURN  is  called  to  calculate  the  burning  rate  perturbations  at  each  of  the  points 
called  for.  Burning  rate  solutions  are  obtained  only  for  every  NTB  th  point  on 
the  initial  line.  The  burning  rate  derivatives  are  then  calculated  by  simple  first 
order  differencing  and  the  mass  burning  rate  at  each  location  is  found  by  adding 
the  perturbation  to  the  initial  value.  The  burning  rate  at  those  points  where  it  is 
not  directly  calculated,  is  then  found  by  linear  interpolation. 

At  this  point  the  solution  is  then  printed  out  every  NPRNT  characteristics 
cycles.  The  pointers  K  and  L,  used  to  direct  storage  for  subroutines  FPT,LHB 
and  RHB,  are  then  shifted  back  to  their  original  values  in  preparation  for  the  next 
cycle. 


Before  the  next  cycle  is  begun,  however,  several  tests  and/or  operations  are 
performed.  If  the  last  calculation  carried  the  solution  past  the  next  plot  time,  the 
variables  of  interest  are  stored  for  plotting.  If  the  number  of  times  this  information 
has  been  stored  (JPL0T)  is  equal  to  the  maximum  specified  (NPLTF),  subroutine 
D0PL0T  is  called  and  the  stored  information  is  plotted  (NPLTF  curves  to  each  graph). 

In  order  to  reduce  the  storage  required  to  compute  the  burning  response  over 
many  cycles  part  of  the  past  history  is  dropped,  periodically.  After  NINT  bum  rate 
calculations  have  been  performed  the  burning  rate  perturbation  is  calculated  as  the 
sum  of  two  parts  (see  subroutine  TRBURN) .  The  first  part  due  to  a  disturbance 
starting  at  t  =  o  and  ending  at  the  time  corresponding  to  the  NINT  th  calculation 
(*NINT)  the  second  part  due  to  a  pressure  and/or  velocity  disturbance  which  originates 
at  tNINT.  After  2  *  NINT  calculations  the  effect  of  the  disturbances  during  the 
first  interval,  up  to  *NINT,  is  dropped  completely  (it  should  have  become  negligible 
if  NINT  has  been  chosen  large  enough).  This  process  Is  accomplished  by  storing 
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the  second  interval  disturbances  (PTB2  and  VTB2)  and  their  corresponding  bum 
rate  perturbations  (EMTB2  and  EMTBV2)  over  their  respective  first  interval  values. 
The  storage  for  the  time  vector,  TIMTB  and  time  interval  array,  DTIM,  is  also 
shifted,  and  the  counter  NCYC  is  reset  to  NINT  (it  was  2  *  NINT).  The  values  of 
the  running  integrals  (RUNI  etc.)  are  also  reset.  This  process  is  then  replaced 
every  additional  NINT  cycles.  Hence,  if  desired,  large  times  can  be  reached 
without  running  out  of  storage. 

Before  the  storage  shifts  are  effected  pertinent  variables  from  the  first  in¬ 
terval  are  printed  out  before  the  information  is  lost.  (Future modification  could 
allow  this  information  to  be  saved  on  tape,  disk  or  drum  files). 

After  the  above  portion  of  the  routine  the  current  time  is  compared  to  TMAX. 

If  TMAX  has  not  been  reached  the  whole  integration  procedure  is  repeated  r.gain. 

If  TMAa  has  been  reached  the  program  checks  to  see  if  any  plots  versus  axial  dis¬ 
tance  have  been  stored  and  not  plotted.  If  so,  D0PL0T  is  called.  Additionally,  if 
lr  ->n  NINT  burn  rate  calculations  have  been  made  when  TMAX  is  reached  PTB1 

o., '  . fBl  are  plotted  versus  TIMTB  atx  =  o  and  x  =  1.  If  velocity  coupling  is 

being  considered  VTB1  and  EMTBV1  are  then  plotted  versus  TIMTB  at  several  dif¬ 
ferent  x  locations. 


SUBROUTINE  INTERP 

After  each  cycle  of  the  characteristics  solution  is  completed  this  subroutine 
interpolates  within  the  characteristics  mesh  to  obtain  a  set  of  points  at  the  original 
x  locations  on  the  initial  line,  all  at  the  same  time  (TMIN). 

As  a  result  of  the  mesh  geometry  and  the  velocity  at  the  end  of  the  grain,  the 
right  hand  boundary  point  always  occurs  at  a  time  somewhat  less  that  all  of  the 
other  characteristics  intersections.  The  routine,  therefore,  extrapolates  (linearly) 
to  establish  a  right  hand  boundary  values  at  TMIN. 


I 

Linear  interpolation  is  then  used  to  find  values  for  the  variables  of  in-  j 

l 

terest  at  each  of  the  original  x  locations.  At  this  stage  each  of  the  points  will,  j 

in  general,  be  at  a  different  time.  Another  linear  interpolation  is  then  perform¬ 
ed  to  obtain  a  series  of  points,  all  at  TMIN.  j 

I 


SUBROUTINE  ITER  (FI,  XI,  XNEW,  NOO) 

This  subroutine  finds  the  root  of  an  equation  by  the  secant  method. 


SUBROUTINE  LEGS 


This  subroutine  is  part  of  the  set  of  steady  state  routines.  It  is  a  general 
matrix  inversion  routine  and  called  from  subroutine  NEWTI  to  aid  in  calculating 
the  next  guess  for  the  chamber  pressure,  PF. 


SUBROUTINE  LHB 

This  subroutine  is  part  of  the  method  of  characteristics  solution  and  solves 
for  points  on  the  left  hand  boundary  x  =  o.  It  follows  the  analytical  description 
of  the  solution  faithfully,  and  as  such,  does  not  require  detailed  discussion.  (See 
subroutine  FPT  for  some  additional  comments). 
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SUBROUTINE  LINEX 


This  subroutine  performs  additional  computations  required  before  the  char¬ 
acteristics  solution  can  be  initiated,  as  follows: 

1.  If  NSET=  1,  the  steady  state  solution  is  bypassed  and  uniform  initial  con¬ 
ditions  with  no  mean  flow  arc  act  up  in  subroutine  LINEX.  This  option  is  used 
only  under  special  circumatances  when  the  program  is  being  applied  to  something 
other  than  a  solid  rocket  motor. 

2.  The  steady  state  solution,  calculated  previously,  is  modified  by  the  addition  of 
a  pressure  and/or  a  velocity  perturbation.  Currently,  the  initial  perturbations  are 
of  the  form 

P  =  Pc.  +  AP  cosfcx) 
v  =  v  +  ..'V  sin(nx) 

P  and  P  are  normalized  by  the  chamber  pressuie,  PF;  while  v  and  Ay  are  normalized 
by  the  chamber  sound  speed,  A.P.  In  the  Luturo,  the  orogram  should  be  modified  to 
3’1-y.v  a  choice  between  severaj  types  of  initial  perturbations;  or  the  perturbation 
calculation  could  be  made  into  a  Program  function,  so  as  to  bo  easily  modifiable. 

The  initial  gas  temperature  and  density  are  modified  i  sen  tropically,  to  reflect 
the  initial  pressure  change.  The  initial  particle  /elocity,  temperature  and  density 
are  currently  not  perturbed. 


3.  The  next  portion  of  the  subroutine  performs  a  series  of  initial  calculations  re¬ 
quired  for  the  transient  burning  rate  calculation.  First,  blank  common  is  zeroed. 
The  transient  burning  rate  parameters  A,  B  and  Bv  are  then  tested  to  see  if  they 
are  compatible  with  the  constraint 


B  +  1 

Any 


If  not,  the  solution  is  terminated.  The  remainder  of  the  routine  either  cal¬ 
culates  initial  values  of  quantities  at  t  =  o,  or  quantities  which  need  only  be  cal¬ 
culated  once. 
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PROGRAM  MAIN 


The  main  program  calls: 

Subroutine: '  Input  to  read  in  and  process  the  input  data; 

Subroutine:  STEDYST  to  solve  for  the  steady  state  solution; 
Subroutine:  LINEX  to  prepare  the  initial  line  data  and. 

Subroutine:  INT  to  carry  out  the  method  of  characteristic  solution. 


SUBROUTINE  MPL0TS 


This  subroutine  is  a  general  CALCOMP  plot  routine  and  is  used  as  follows: 
CALL  MPL0TS  (X,  Y,  N,M,N0,LEN,  TTTLX,  TITLY,  HEAD,  NX,NY,NH,FAC,LABL) 
where 

X  =  vector  or  doubly  dimensioned  array  (N  x  1  minimum)  of 

independent  variables  (see  M  for  usage) 

Y  =  doubly  dimensioned  array  of  the  dependent  variables  to  be 

plotted.  (N  x  N0  minimum) 

N  =  actual  dimensioned  length  of  the  columns  for  the  X  and  Y 

arrays.  N  must  be  at  least  2  cells  longer  than  the  used 
dimensions  of  the  X  and  Y  arrays. 

M  =  a  flag  such  that  if 


M  >0 

there  is  an  X  vector  corresponding  to  each 
Y  vector  curve,  i.e.,  the  Y(1,J)  vector 
would  be  plotted  versus  the  X(1,J)  vector. 

M<0 

there  is  only  one  X  vector  corresponding  to 
all  of  the  Y(1,J)  curves. 

M  =  0 

same  as  M  >0  except  no  legends  will  be 
plotted  out. 

1 M  1  =  1 

legends  will  be  plotted  out  with 

LABL(l)  =  Hollerith  header 

LABL(J)  =  numerical  value  corresponding 

to  the  J-l  curve 

!  M  |  =  2 

legends  will  be  plotted  out  with 

LABL(I)  =  Hollerith  name  for  the  Ith  curve 
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N0 

LEN(I) 


TITLX(I)  = 
TITLY(I)  = 
HEAD  (I)  = 


NX 

NY 

NH 

FAC 


I,ABL(I) 

Common  Linkage 


•  M  |=  3  no  legends  will  be  plotted  out 

.  M  1=  10,20  or  30  same  as  ;_Mj  =  1,2,  or  3  except  that 

a  border  will  be  put  on  the  plot  frame 

number  of  curves  to  be  plotted  per  frame 

a  vector  of  length  NO  containing  the  length  of  the  Ith 
curve  (i.e.,  number  of  data  points  stored  in  the  Ith 
column  of  the  Y(J,I)  array). 

a  vector  containing  the  Hollerith  label  for  the  X-axis. 

a  vector  containing  the  Hollerith  label  for  the  Y-axis. 

a  vector  containing  the  Hollerith  label  for  the  master 
hooding. 

number  of  characters  in  the  X-axis  label. 

number  of  characters  in  the  Y-axis  label. 

number  of  characters  in  the  master  heading  label. 

the  factor  by  which  the  plot  will  be  scaled.  FAC  =  .8 
will  pioduce  a  plot  compatible  with  an  8  1/2  x  11  page 
size. 

vector  containing  either  Hollerith  and/or  floating  point 
numerical  data  (see  M -description). 


In  addition  to  the  calling  sequence,  a  labeled  common  area  exists  through 
which  the  user  can  exercise  additional  control  over  the  plots.  The  label  common 
area  is 

C0MM0N/RR1L/AXL ,  AYL ,  FT,  ISYM ,  INC ,  LNTYP ,  NLPC , NDEC 

where 


AXL  =  10.0  lenght  of  the  X-axis 

AYL  =  8.0  lenght  of  the  Y-axis 

HT  =  .  14  character  height  for  header  labels 

ISYM  =  4  1SYM+1  is  the  number  for  the  starting  symbol  on  the  graphs 

INC  =  1_  increment  at  which  data  points  within  a  vector  will  be  plotted. 

LNTYP  =  _1  see  the  CALC0MP  write  up  for  subroutine  LINE  for  the  usage 

of  this  variable 

NLPC=  4  number  of  curve  legends  per  column 

NDED  =  5  used  v>  hon  !  M  !=  1  or  10.  Controls  the  number  of  digits  to 

the  right  of  <.he  decimal  point  on  the  values  usee  in  the  legends. 

Note:  The  values  underlined  are  the  nominal  values  and  have  boon  set  via  DATA 

STATEMENTS.  Therefore  in  order  to  change  those  values,  you  should  use  an  execute 

statement  in  your  program. 
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SUBROUTINE  NEW'fl 

'iliis  subroutine  is  part  of  the  next  of  subroutines  for  solving  the  steady 
state  two  point  boundary  value  problem.  It  uses  Newton's  method  to  generate 
new  guesses  for  those  initial  values  which  must  established  in  order  to  satisfy 
the  boundary  condition(s).  In  the  current  application,  NEWT1  generates  guesses 
for  the  chamber  pressure,  PF,  until  the  mach  number  at  the  end  of  the  grain  matches 
that  calculated  from  a  fractional  lag  nozzle  solution. 


SUBROUTINE  PHIF 

This  subroutine  is  called  from  subroutine  FRAKLAG  and  is  used  in  conjunction 
with  subroutine  ITER  to  insure  that  the  calculated  values  of  K  and  GBAR  (y)  are  com¬ 
patible. 

If  y  =  y  there  are  no  particles,  and  the  flag  KTEST  is  set  equal  to  1.  This 
flag  is  used  elsewhere  in  the  program  to  determine  which  calculations  may  be  skip¬ 
ped  when  there  are  no  particles . 


SUBROUTINE  PRNT 

This  subroutine  prints  out  the  initial  line,  as  well  as  the  solution. 

0  After  each  NPRNT  cycles  of  the  characteristics  calculation. 

0  At  each  time  when  the  variables  are  stored  for  plotting. 

0  At  the  final  time. 

The  header  card  containing  case  information  is  printed  at  the  top  of  each  page. 
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FUNCTION  RATE  (PF  .RI,  UI .  TI,G,  A) 


This  function  computes  the  steady  state  mass  burning  rate,  ,  in  nondimen- 
sional  form. 


SUBROUTINE  RHB 


This  subroutine  is  part  of  the  method  of  characteristics  solution  and  solves 
for  p  lints  on  the  right  hand  boundary  x  =  1.  It  follows  the  analytical  descr  otion 
of  the  solution  faithfully,  and  does  not  require  detailed  discussion.  (See  sub¬ 
routine  FPT  for  some  additional  comments). 


FUNCTION  RHS  (FG,  A.RHO,  RHOP.  U,  UP .  T,  TP  .W.DTIM  ,X) 


This  function  is  called  by  FPT,  RHB  and  LHB  to  evaluate  the  right  hand  sides 
of  the  compatibility  relations  for  both  the  momentum  and  energy  equations. 

The  parameter  FG  j?  used  to  control  which  characteristics  line  is  being  con¬ 
sidered,  as  shown  below: 

FG  =  -1  Right  running  characteristic 

FG  =  1  Left  running  characteristic 

FG  =  2  Energy  equation  along  gas  streamline 
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SUBROUTINE  RKAM 


This  subroutine  is  part  of  the  set  of  subroutines  which  solves  the  steady 
state  two  point  boundary  value  problem.  RKAM  is  a  general  subroutine  for  in¬ 
tegrating  simultaneous  ordinary  differential  equations.  It  allows  a  choice  be¬ 
tween  the  Adams-Moulton  and  4th  Order  Runge-Kutta  methods.  In  the  present 
program  the  Adams-Moulton  method  is  used  to  integrate  the  six  simultaneous 
differential  equations  which  describe  the  steady  state  flow  In  the  chamber. 


SUBROUTINE  SCAL 


Subroutine  SCAL  is  used  in  conjunction  with  the  plot  subroutine  MPLOTS. 
It  is  used  to  optimize  the  selection  of  a  plot  scale . 


SUBROUTINE  3TEDYST 


This  subroutine  does  the  following: 

1.  Call  subroutine  FRAKLAGto  find  the  fractional  lag  nozzle  solution. 

2.  Solves  a  transcedental  equation  to  find  the  mach  number  at  the  end  of  the  grain 
as  a  function  of  port  to  throat  area  ratio.  (Uses  the  results  of  the  fractional 
lag  solution). 

3.  Sets  the  boundary  values  at  x  =  o  for  each  of  the  six  variables  FI(1)  . . .  FI(6). 

4.  Calculates  chamber  fractional  lag  parameters,  AK  and  AL.  These  parameters 
are  used  to  obtain  an  approximate  solution  for  the  particle  flow  in  the  chamber  when 
the  flow  is  near  dynamic  equllibirum,  and  hence  is  subject  to  numerical  difficulties 
when  treated  in  the  regular  manner. 
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5.  Subroutine  BVP  (which  is  actually  a  nest  of  subroutines)  is  called  to  solve 
for  the  steady  state  flow  conditions  in  the  chamber. 


SUBROUTINE  TFUNC 


This  subroutine  computes  gas  viscosity  as  a  function  of  temperature.  The 
relationship  currently  contained  in  the  program  is  actually  Sutherland's  law  for 
air.  This  functional  relationship  should  be  changed  to  more  accurately  reflect  the 
viscosity  of  the  gas  mixtures  found  in  solid  rocket  motors. 


SUBROUTINE  TRBRNA(TMIN) 

Each  time  a  new  set  of  transient  burning  rate  calculations  is  to  be  made  this 
subroutine:  updates  the  required  counters;  stores  the  current  time  in  the  TIMTB 
array;  calculates  the  size  of  the  latest  interval,  DELTA,  and  stores  it  in  the  DTIM 
array.  The  TEM(I)  array  used  in  the  trapezoidal  integration  procedure,  is  also  cal¬ 
culated. 


SUBROUTINE  TRBURN(IK) 

This  subroutine  calculates  the  transient  burning  rate  response  at  the  IK  th  x 


location. 


The  subroutine  is  actually  divided  into  two  distinct  parts.  Each  part  is 
structured  identically:  the  first  part  calculates  pressure  coupled  response;  if 
KPRES  /  0;  the  second  part  calculates  velocity  coupled  response,  if  KVEL  /  0. 

Each  of  these  parts  is  further  divided  into  subsections.  When  MCYC  < 
NINT  the  burning  rate  is  calculated  normally  as  a  function  of  a  single  pressure 
(PTB1)  and/or  velocity  (VTB1)  perturbation. 


MCYC  <  NINT 

{only  pressure  coupling  is  illustrated) 


EMTB1 


NINT 

Number  of  cycles 


2NINT 


NINT 

Number  of  cycles 


2NINT 


When  MCYC  is  between  NINT  and  2  NINT  the  burning  rate  is  calculated  as 
the  sum  of  the  response  to  two  perturbations,  as  shown  below. 
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MCYC  >  NINT 


{only  pressure  coupling  is  illustrated) 


EMTB1 


NINT 


2NINT 


NINT 


2NINT 


EMTB2 


NINT 


2NINT 


NINT 


2  NINT 


Number  of  cycles 


Number  of  cycles 


In  an  effort  to  minimize  the  computation  time  associated  with  the  calcu¬ 
lation  of  the  burning  rate,  the  numerical  evaluation  of  the  quadratures*  is  broken 
into  two  parts .  The  integrals  in  which  t  is  a  parameter  must  be  evaluated  from 
t  »  0  to  t  =  tm^n  -  6,  each  time  (unless  part  of  the  past  history  has  been  dropped 
(see  subroutine  INT).  These  integrals  are  calculated  using  the  Trapezoidal  Rule. 
The  integrals  involving  £  only  can  be  evaluated  In  a  running  manner;  the  contribu¬ 
tion  from  each  successive  interval  being  added  to  the  sum  of  the  previous  ones 
(RUNI,  RUN VI  etc.). 


FUNCTION  VPERT  (ARG1.ARG2) 

The  calculation  of  the  velocity  coupled  transient  bum  rate  perturbation  re¬ 
quired  the  following  expression  to  be  evaluated 

i^lul  -  ut)  -  v2(u  -  ut) 

WithARGl  =  lul  and  ARG2  =  u{  the  first  part  of  this  expression  is  evaluated. 
ARG1=  Uo  and  ARG2  =  ut  yields  the  latter  part  of  the  expression. 


See  Equation  (4-68)  in  volume  1. 
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3.  DICTIONARY  OF  COMMON  VARIABLES 

All  of  the  variables  which  appear  in  Common*  (blank  or  labeled! 
are  defined  in  this  section.  Blank  Common  is  considered  first,  followed  by 
the  labeled  Commons  in  alphabetical  order. 


COMMON 

BLANK 


NAME 


PTBl(IJ) 


PTB2(I,  J) 


EMTBl(IJ) 

EMTB2(I,J) 

TIMTB(I) 

DTIM(I) 

TEM(I) 

EMTB(I.J) 


MATH.  SYM, 

(£) 

F  2 

( —  ) 


DXTB 

TCtfN(I) 

VTBl(IJ) 

VTB2(I,I) 

EMTBVl(IJ) 

EMTBV2  (I,  J) 


r/4x 


DESCRIPTION 
first  pressure  perturbation 

second  pressure  perturbation 

burning  rate  response  to  first 
pressure  perturbation 

burning  rate  response  to 
second  pressure  perturbation 

time  of  Ith  set  of  burning  rate 
calculations 

TIMTB (I)  -  TIMTB(I-l) 

multiplier  used  in  calculating 
burning  rate  integrals 

total  burning  rate  perturbation , 
J=1  next  to  last  time,  1=2 
current  time 

total  number  of  X  locations 
at  which  transient  burn  rate  is 
calculated 

distance  between  locations  at 
which  transient  bum  rate  is 
calculated 


[c i  (|u|  -  Ut)  -C  2  (U  -  Ut)]i  first  velocity 

perturbation 

[ei  (|  U|  -  Ut)  -c2  (U  -  Ut)]9  second  velocity 

perturbation 


m  o 


burning  rate  response  to  first 
velocity  perturbation 

burning  rate  response  to 
second  velocity  perturbation 


♦With  the  exception  of  some  common  blocks  used  entirely  within  the  confines  of  the 
BVP  rest  of  subroutines  for  solving  the  steady  state  two  point  boundary  value  problem. 


*73W*W*fP 


COMMON 

BLANK 


ACE 


BNDS 


NAME 


MATH  SYM , 


DESCRIPTION 


CI(I),  C2(I) , 

03(1),  C4(J) 

constants  in  pressure  coupled 
burning  rate  integral 

CV1(I),  CV2(I) , 

CV3(I),  CV4(I) 

constants  in  velocity  coupled 
burning  rate  integral 

RUN  1(1) 

part  of  pressure  coupled  burning 
rate  response  integral  that  is 
evaluated  in  running  manner 

RUN  2(1) 

see  RUN1  but  for  second 
pressure  disturbance 

RUNVl(I) 

part  of  velocity  coupled  burning 
rate  response  integral  that  is 
evaluated  in  running  manner 

RUNV2  (1) 

see  RUNV1  but  for  second 
velocity  disturbance 

DWDT(I)  da1 

dt 

time  derivative  of  burning 
rate 

This  common  contains  all  of  the  quantities  which  can  be  input 
to  the  steady  state  solution  package  (subroutine  BVP).  Those 

of  interest  in  the  current  program 

of  input. 

are  discussed  in  the  description 

K 

storage  pointer  in 
characteristics  solution 

L,  ELL 

storage  pointer  in 
characteristics  solution 

NS1 

set  equal  to  1 

NE1 

set  equal  to  N 

NS2 

set  equal  to  2 

NE2 

set  equal  to  N 

FLAG 

not  used 

N 

M-l 

M 

number  of  axial  locations 

in  characteristics  mesh 


I 


I 
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E 

1 


i 


S 

c 


COMMON 

CHAR 


FLOS 


FPL0T 


tmtiaSaiUtialU 


name  math  sym. 

TIM(IJ)  t 


x(i,J) 

X 

U(I,J) 

U 

UP(IJ) 

Up 

T(I,  J) 

T 

TP(IJ) 

Tp 

RH0(IJ) 

0 

RHOP(IJ) 

°p 

p(I/ J) 

p 

W(IJ) 

60 

WP(I,J) 

% 

A(I,J,ADUM(IJ) 

a 

UI  u 

rap  Up 

RH0I  0 

RH0IP  Qp 

U  T 

TIP  Tp 


NFLTF 

TPL0T(I) 

HEADER(I) 
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DESCRIPTION 

time 

axial  distance 
velocity  point 
particle  velocity 
temperature 

particle  temperature 

*  / 

density 

particle  density 
pressure 

mass  burning  rate 
particle  mass  burning  rate 
sound  speed 


velocity 

particle  velocity 
density 

particle  density 
temperature 
particle  temperature 

maximum  number  of  plots  that 
can  be  drawn  on  one  curve 

plot  times  (see  input  descrip, ) 

contains  case  title  (see  input 
de  scrip.) 


a*  ^ -  - 


COMMON  NAME 


FRAK 


ME8 


MATH  SYM. 


DESCRIPTinM 

mach  number  of  equivalent  gas 


relates  M  to  M 
e  e 


FSTART 


FI(6) 


ifonol  °f  ^teady  stQte  variables 
(conserved  quantities) 

chamber  pressure  (psf) 


MCYC 


NCYC 


NKM2 

NCYC2 


NK2M2 


DELTA 


SQDEL 


NPRNt,  uth,  NINT,  KVEL , 

kpres,cnv,bv,dpv 


number  of  times  transient  burn¬ 
ing  rate  has  been  calculated 

MCYC  +  1 

number  of  cycles  over  which 
burning  rate  integral  must  be 
integrated 

NCYC  +  1 

NK  -  2 

NCYC  -  NINT 
NCYC2  +  1 
NK2  -  2 

indicates  which  x  location  is 
being  considered 

time  interval  between  current 
and  last  burning  rate  solution 


see  description  of  input 


’  n  int  iwi  »>■'.*  yij*!1;?  y  ■'  m  .y 


!i 


i^1 


5S5* JW“ ^V.l t* L  .' '! .«reW# ^ 


common 

NTERP 


PL0TQ 


NAME 

MATH  SYM. 

DESCRIPTION 

All  of  these  arrays  contain  the  quantities  calculated  at  the  end  1 

of  the  second  characteristics  sweep,  but  before  Interpolation  j 

TIMT(I) 

time 

XI{I) 

axial  distance 

IJT(I) 

velocity 

UPT(I) 

particle  velocity 

TT(I) 

temperature 

TPT(I) 

particle  temperature 

RH0PT(I) 

density 

PT(I) 

pressure 

WT(I) 

mass  burning  rate 

TMIN 

time  at  the  end  of  the  last 
characteristics  cycle 

XPL0T(I,J) 

axial  distance 

PPL0T(I,  J) 

p-F 

pressure  perturbation 

UPL0T(I,J) 

u-u 

velocity  perturbation 

UPPL0T(IJ) 

Up -Up 

particle  velocity  perturbation 

WP10T(I,J) 

w-W 

mass  burning  rate  perturbation 

LEN(I) 

number  of  points  stored  for 
plotting 

XLABL(I) 

see  subroutine  MPL0TS 

NH 

number  of  Hollerith  characters 
in  Header  Title  j 

JPL0T 

\ 

number  of  curves  stored  for  plot-  • 
ting  on  each  graph 

1  -i 


COMMON 


NAME 


MATH  SYM, 


DESCRIP  HON 


PR0X 


velocity  fractional  lag 
parameter  in  chamber 


temperature  fractional  lag 
parameter  in  chamber 


chamber  pressure  (psf) 


RR1L 


see  subroutine  MPL0TS  description 


TESTING 


KTEST 


KTEST  =  1  if  there  are  no  par¬ 
ticles 


EPS1 , EPS2 


*l'e2 


convergence  criteria  for  char¬ 
acteristics  solution 


TZER0 


PZER0(I) 

UZER0(I) 

UPZER0(I) 


f  +  pi  =  o 

u+uUo 


perturbed  pressure  at  t  =  0 
perturbed  velocity  at  t  =  0 
particle  velocity  at  t  =  0 


WZER0(I) 


burning  rate  at  t  =  0 


COMMON 

VARS 

ZJN 


MATH  SYM.  DESCRIPTION 

see  Description  of  Input 

IREAD  see  Description  of  Input 


4.  DESCRIPTION  OF  PROGRAM  INPUT 

Program  input  for  the  instability  program  consists  of  two  distinct 
groups.  The  first  set  of  input  controls  the  instability  solution,  while  the 
second  set  controls  the  steady  state  solution  which  provides  initial  conditions 
for  the  instability  solution.  Many  of  the  input  quantities  have  been  assigned 
preset  nominal  values;  these  quantities  need  be  input  only  if  it  is  desired  to 
alter  the  preset  value.  Nominal  values  are  indicated  where  applicable.  In 
many  cases,  the  second  set  of  input  can  be  eliminated  completely,  since 

nominal  values  have  been  preset  for  all  of  these  quantities  (see  discussion 
of  IREAD,  below). 


4 . 1  Instability  Program  Input 

The  first  card  of  this  set  contains  the  case  title  in  columns  1-70.  This 
title  will  appear  at  the  top  of  each  page  of  output  and  on  each  computer  plot. 
The  input,  itself,  follows  the  title  card,  and  follows  the  standard  NAMELIST 
format.  The  first  card  of  this  set  must  contain  $DATA  starting  in  column  2. 
The  last  card  must  contain  $END  starting  in  column  2. 

The  input  variables  in  the  DATA  namelist  are  described  below  in 
functional  sub-groups. 


Engine  Geometry 
DP0RT 
DTHR0T 
L 

RC 

SK 


Chamber  port  diameter,  Dp,  in  inches. 

Nozzle  throat  diameter,  Dt<  in  inches. 

Length  of  the  grain,  L,  in  inches. 

Nondimensional  throat  radius  of  curvature,  R, 
Rc/rf  (Nominal  value  =1.0). 

Area  =  SK(X-1)  +  Ae,  currently  only  SK=0 
(constant  area)  should  be  used. 
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Steady  State  Burning  Rata  Cnna^nt. 

The  steady  state  burning  rate  is  specified 
CTILDA 
CK 
CN 
PREF 


as  r  =  C(P/Pref)n(i  +cku) 

Constant  in  burning  rate  expression,  C,  in/sec. 
Erosive  burning  constant,  C^,  sec/ft. 

Pressure  exponent,  n. 

Normalizing  pressure,  psia  (nominal  value  =  500) 


Gas  Constanta 
CP 

Gamma 

Pr 

TF 

Propellant  Constant 
RH0S 
RH0M 

CKAPA 


Gas  specific  heat,  Cp,  Btu/ib  °R. 

Gas  isentropic  exponent,  ,  . 

Prandtl  number,  Pr,  (Nominal  value  =  1.0) 
Adiabatic  flame  temperature,  Tp,  °R 


Density  of  the  solid  propellant,  t  ,  g/cc. 

s 

Density  of  the  metal  oxide,  ,  g/cc 
(Nominal  value  =  4),  m 

Thermal  diffusivity  of  the  propellant, 

Rs ,  cnrvsec. 


i 


PDIA 

A 

B 


C 

CNV 


t 


Particle  diameter  in  microns. 

Transient  burning  rate  parameter. 

Transient  burning  rate  parameter. 

Velocity  coupling  transient  burn  rate 
parameter  (currently  set  equal  to  B  in 
Subroutine  Input) . 

Ratio  of  particle  to  cas  specific  heats, 
Cs/CP. 


Velocity  coupling  transient  burn  rate 
parameter  (currently  set  equal  to  CN 
in  Subroutine  Input). 


BETA1 


UTH 


Program  Constants 
H 

IM 

EPS1 

EPS  2 

TMAX 


PFO 

DP 

DPV 


Flags  and  Counters 
NPRNT 

MAXIT 


SHmiast 


Particle  to  gas  weight  flow  ratio,  B, . 

(To  run  a  case  with  no  particles  requires 
only  that  B^  be  set  equal  to  zero) . 

Threshold  velocity  to  erosive  burning 
response  (Nominal value  =0). 


Integration  step  size  for  steady  state 
solution  (Nominal  value  =  .01). 

Total  number  of  points  at  which  steady 
solution  is  obtained,  =  1/H+l  (Nominal 
value  =  101) . 

Convergence  criteria  for  p,  and  a,,  c, . 
(Nominal  value  =  1  x  lO-5);  6  1 

Convergence  criteria  for  Op  *  €«?•  (Nominal 
value  =  1  x  10-5).  *3  c 

Final  value  of  time.  Computation  ceases 
when  t  >T  .  (Nominal  value  =  0.1) . 
Currently  tnex value  of  TMAX  must  be  selected 
such  that  not  more  than  500  of  the  character¬ 
istics  calculation  cycles  (not  wave  cycles) 
will  be  computed.  In  the  nondimensional 
coordinate  system  used ,  the  time  increment 
for  one  computation  cycle  is  approximately 
equal  to  the  distance  between  points  on  the 
initial  line  -  For  this  reason  the  total  compu¬ 
tation  time  goes  up  with  the  square  of  the 
number  of  points  on  the  initial  line . 

Initial  guess  for  chamber  pressure,  psia. 

Magnitude  of  the  initial  pressure  pulse, 
AP/PF. 

Magnitude  of  the  initial  velocity  pulse, 
AU/aj.  (Nominal  value  =  0) . 


Solution  is  printed  every  NPRNT-th  calcula¬ 
tions!  cycle ,  (Nominal  value  =  10) . 

Maximum  number  of  iterations  allowed  for 
the  characteristics  solution  at  a  given  point. 
The  solution  continues ,  even  if  it  has  not 
converged.  •  lominal  value  =  3). 
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J 

V 

\ 

f> 

y 


J 


l tur'WT-ti 


Provides  an  option  for  bypassing  the  steady 
state  solution.  NSET  =  0  is  the  nominal 
value,  and  a  steady  state  solution  is  obtained. 
NSET  =  1  provides  for  the  generation  of  uniform 
initial  conditions .  This  option  should  not  be 
used.  It  was  set  up  specifically  for  some  check 
out  runs  and  is  not  now  currently  operational. 

This  flag  determines  whether  the  $INPUT 
NAMELIST  SET  must  be  included  in  the  input 
deck.  If  IREAD  =  1  (Nominal  value),  the 
$INPUT  Namelist  data  is  not  required  and 
no  attempt  will  be  made  to  read  it.  The 
nominal  value  specified  for  each  of  the 
quantities  is  used  by  the  program.  Use 
IREAD  =  0  ,  if  it  is  desired  to  change  one 
or  more  of  the  $INPUT  nominal  values .  With 
IREAD  =  0,  the  $INPUT  is  read,  in  Subroutine 
BVP. 

Sets  the  maximum  number  of  curves  that  will 
be  drawn  on  one  computer  plot.  (Nominal 
value  =  4,  can  be  set  less  than  4,  but  not 
higher,  unless  the  necessary  array  dimensions 
are  increased). 

The  steady  state  solution  at  every  NTRIG-th 
point  is  stored  for  use  on  the  initial  line  for 
the  method  of  characteristics.  The  maximum 
total  number  of  points  allowed  on  the  initial 
line  is  200.  The  number  of  points  on  the 
initial  line  is  given  by  M  =  1+IM/NTRIG, 
done  in  integer  arithmetic.  NTRIG  has  been 
assigned  a  nominal  value  of  2,  which  provides 
51  points  on  the  initial  line  when  IM  =  101. 

The  transient  burning  rate  integrals  are  evalu¬ 
ated  starting  at  every  NTB-th  location  on  the 
initial  line.  Burning  rate  values  at  other 
locations  are  computed  by  interpolation.  The 
total  number  of  points  on  the  initial  line,  less 
one,  must  be  exactly  divisible  by  NTB,  i.e. , 
(M-1)/NTB  must  be  a  whole  number.  The 
program  will  automatically  terminate  if  there 
is  an  improper  correspondence  between  the 
values  of  M  and  NTB.  The  maximum  number 
of  locations  at  which  the  burning  rate  integrals 
may  be  evaluated  is  currently  26. 

After  one  transient  burning  rate  calaulations 
have  been  performed  NINT  times ,  the  burning 
rate  is  calculated  as  the  sum  of  one  response 
to  two  perturbations.  After  2  NINT  calculations 
the  response  to  the  first  disturbance  is  deleted 
and  storage  is  shifted.  (Nominal  value  = 
1,000,000). 


KVEL 


fl)  velocity  coupling  included ,  =  0  no 
velocity  coupling  (Nominal  value  =  0) . 


KPRES  ^0  pressure  coupling  included,  =  0  no 

pressure  coupling  (Nominal  value  =  1). 

In  addition  to  the  above  quantities ,  an  array  must  be  specified  to  govern 

the  times  at  which  the  variables  will  be  stored  for  computing  plotting . 

This  array  is  denoted  by  TPL0T(I)  and  has  the  maximum  dimension  of  20. 

The  variables  P'  =  P-P  ,  u'  =  u-u  ,  u'  =  u  -u  and  <*>*  =  Ok  are  stored 

o  o  p  p  pQ 

at  each  x  location,  the  first  time  a  calculation  cycle  ends  at  a  time 
greater  than  or  equal  to  each  successive  value  of  TPL0T(I).  Each  time 
the  variables  have  been  stored  at  NPLTF  different  times  Subroutine  MPL0TS 
is  called;  which  then  processes  the  stored  information  and  writes  out  the 
plotting  instructions  on  tape.  If  19  or  less  plot  times  are  specified  the 
program  automatically  inserts  t  =  1,000,000  as  the  last  plot  time.  There¬ 
fore,  if  plotting  is  not  desired  TPL0T  need  not  be  input,  as  TPL0T(1)  will 
be  set  to  ..,000,000  and  will  never  be  reached. 

If  TPL0T(1)  £  99  the  program  also  plots  the  pressure  and  transient 
burning  rate  histories  at  x  =  0  and  x  =  1 .  One  can  choose  to  have  these 
latter  plots,  and  not  the  aforementioned  set,  by  choosing  a  value  for  TPLOT(l) 
greater  than  TMAX,  but  less  than  99. 

Steady  State  Input 

The  quantities  which  control  the  steady  state  solution  are  also  read  in 
using  the  NAMELIST  format.  The  first  card  of  this  set  should  read  $INPUT, 
and  the  last  card  $END.  When  used,  the  $INPUT  data  immediately  follows 
the  $END  card  of  the  $DATA  NAMELIST.  The  quantities  that  may  be  input  are 
described  below.  Nominal  values  have  been  preset  for  each  quantity  and  are 
indicated. 

HI  Integration  step  size  for  the  steady  state 

solution  (Nominal  value  =  .01). 

LTRIG  Controls  the  locations  at  which  the  steady 

state  integrated  values  and  derivatives  are 
printed.  (LTRIG  controls  the  printing  of  the 
conserved  quantities  and  their  derivatives, 
not  the  flow  variables,  P,  u,  P,  etc., 
themselves) . 

LTRIG  =  3,  (Nominal  value)  print  only  at 
x  =  0  and  x  =  1 . 

LTRIG  =  2,  print  at  every  NTRIG-th  step. 


ITRIG 


Determines  the  frequency  at  which  the  conserved 
quantities,  and  their  derivatives ,  are  printed. 
ITRIG  =  1,  the  results  are  printed  each  iteration. 
ITRIG  =  2 ,  the  results  are  printed  only  after 
the  final  iteration ; 

ITRIG  =  3,  (Nominal  value),  the  results  are 
not  printed. 

NTRIG  Must  be  identical  to  NTRIG  in  the  $DATA 

input  set.  Here  it  governs  the  printout  as 
shown,  for  LTRIG  =  2.  (Nominal  value  =  2). 

MAXIT  Not  the  same  as  MAXIT  in  $DATA  input.  Here 

MAXIT  sets  the  maximum  number  of  iterating 
allowed  for  the  steady  state  solution.  If  the 
Mach  number  at  the  end  of  the  grain  is  not 
matched,  to  within  the  convergence  tolerance, 
in  MAXIT  iterations,  the  solution  is  terminated. 
(Nominal  value  =  10). 

KB  If  KB  =  1  a  bounding  procedure  is  used  in 

conjunction  with  the  iterative  solution.  Other¬ 
wise  KB  =  0  (Nominal  value).  Unless  the  initial 
chamber  pressure  guess,  PFO,  is  off  by  at 
least  an  order  of  magnitude,  or  more,  the 
bounding  option  should  not  be  required. 

BSEVEN  If  KB  =  1,  the  chamber  pressure  will  not  be 

allowed  to  change  by  more  than  the  value  of 
BSEVEN  on  any  one  iteration.  This  bounding 
procedure  will  usually  slow  the  convergence 
rate,  and  therefore  should  not  always  be  used. 
However,  in  some  cases,  particularly  when  the 
initial  guess  is  qui*e  poor,  this  bounding  allows 
a  converged  solution  to  be  achieved  when  the 
iteration  would,  otherwise,  be  unstable. 

A  subset  of  the  subroutines  used  to  solve  the  steady  state  problem 
consists  of  library  subroutines  designed  to  solve  a  general  two  point  boundary 
value  problem  involving  up  to  20  differential  equations  and  constraints.  All  of 
the  options  and  features  of  this  general  routine  are  not  required  in  order  to  solve 
the  present  problem.  As  a  result,  the  list  of  input  has  been  restricted  to  the 
quantities  shown  above.  The  other  items,  normally  required  as  input,  have  all 
been  preset  in  subroutine  BVP  (and  should  not  be  changed),  or  are  of  no 
consequence  in  the  present  program  and  have  been  ignored. 
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5.  DESCRIPTION  OF  PROGRAM  OUTPUT 


The  printed  o.itput  from  the  instability  program  is  described  below.  The 
output  is  presented  in  the  order  in  which  it  appears.  Sample  output  for  the 
test  case  is  contained  in  Section  6  and  should  be  referred  to  in  conjunction  with 
the  text. 

5.1  Printout  of  Program  Input 

The  DATA  namelist  input  quantities  are  printed  out  in  sub-groups  which  cor¬ 
respond  exactly  to  those  shown  in  the  discussion*;  i.e. , 

ENGINE  GEOMETRY 
BURNING  RATE  CONSTANTS 
GAS  CONSTANTS 
PROPELLANT  CONSTANTS 
PROGRAM  CONSTANTS 
FLAGS  COUNTERS 

The  input  units  for  those  quantities  having  dimensions  are  indicated. 

5.2  Fractional  Lag  Solution 

The  next  page  of  output  contains  the  parameters  calculated  as  a  result  of  the 
fractional  lag  nozzle  solution.  The  quantities  themselves  are  easily  '•elated  to 
the  variables  in  the  analytical  description  of  the  solution. 

5.3  Steady  State  Solution 

The  flow  variables  of  interest  are  printed  out  at  the  end  of  each  iteration  of  the 
steady  state  solution.  The  axial  spacing  of  this  printout  corresponds  to  the  spacing 
of  points  on  the  characteristics  Initial  line  and  does  not  represent  all  of  the  points 

calculated,  unless  NTPIG  =  1. 

At  the  top  of  the  page  the  current  value  of  chamber  pressure  is  output  in  psi. 

At  the  end  of  each  iteration  the  number  of  the  Iteration,  and  the  error  in  matching 
the  mach  numbers  at  the  end  of  the  grain,  are  printed  out.  When  convergence  has 

*The  quantities  UTH,  NINT,  KVEL,  KPRES,  CNV,  DPV,  BV  have  yet  to  be  integrated 
into  the  output  format.  These  quantities  are  printed  out  after  the  plot  times. 
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been  achieved  the  message 

STEADY  STATE  SOLUTION  COMPLETED 

is  printed  out. 

5.4  Instability  Solution 

For  each  point  on  the  starting  line  pressure,  velocity,  temperature,  burning 
rate,  particle  velocity,  particle  temperature  and  particle  density  are  printed  out. 
These  quantities  are  also  printed  out  every  time  an  additional  NPRNT  characteris¬ 
tics  cycles  have  been  completed. 

The  case  title,  run  date,  and  time  {in  hours)  are  output  at  the  top  of  each  page, 
to  aid  in  run  identification. 

The  total  number  of  characteristics  cycles  computed  (LCYC)  and  the  corres¬ 
ponding  nondimensional  time,  are  also  printed. 

The  solution  is  also  output  each  time  quantities  are  stored  for  plotting,  as 
well  as  at  the  end  of  the  last  cycle. 

Just  before  the  instability  solution  is  printed  at  the  final  station  a  summary 
of  the  transient  burning  rate  calculations  is  printed.  The  output  quantities  are 


as  follows. 

TIME 

time  at  which  the  results  were  calculated 

PW 

pressure  perturbation 

EMW 

total  pressure  coupled  bum  rate  perturbation 

VW 

velocity  perturbation  (as  used  in  calculating 
velocity  coupled  response) 

EMWV 

total  velocity  coupled  bum  rate  perturnation 

EMW1 

pressure  coupled  response  to  the  first  pressure 
disturbance  only 

EMWV1 

velocity  coupled  response  to  the  first  velocity 
disturbance  only 

! 

i 


i 
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If  KVEL  -  0,  the  aforementioned  quantities  are  printed  only  at  x=0  and  x=l. 
If  KVEL  /  0  transient  burning  rate  results  are  also  printed  at  several  Intermediate 
x  locations. 
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6.  SAMPLE  CASE 


Portions  of  the  output  from  a  sample  case  are  presented  in  this  section  to 
facilitate  program  check-out.  The  sample  case  includes  both  pressure  and 
velocity  coupling. 

Copies  of  the  CALCOMP  plots  generated  by  the  sample  case  are  also  in¬ 
cluded  so  the  operation  of  the  plotting  routines  can  be  checked  out. 
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CONVERGENT  SOLUTION 
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